/*
Copyright (c) 2013 Elizabeth Aston

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/



/* Steady state genetic algorithm
   Input arguments distanceBetweenGenotypes, fitness0, fitness1, radius0, radius1 
   ./GeneticAlgorithm ...arguments... > output.txt
*/


#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

int numGenotypes;
const int genotypeLength = 30;
const int numGenerations = 10000;
const int numMutationsPerReproduction = 1;
int mutationRate; // probability (X in 1000000) of mutation
const int alphabetSize = 4; // Numerical representation of GATC
int **genotypes; /*[numGenotypes][genotypeLength]*/ // Stores the genotypes that make up the population
double *fitness;/*[numGenotypes]*/ // Stores the fitness of each genotype in the population
double fitness0 = 15;
double fitness1 = 10;
int radius0 = 2; //} radius of X means up to X mutations away
int radius1 = 5; //}
int target0[genotypeLength];
int target1[genotypeLength];
int distanceBetweenTargets = radius0+radius1 + 3; // keep in range [radius0+radius1,genotypeLength]
bool peak0;
bool peak1;

const int probabilityOfReplacing1stFittestOf3 = 250; //         1 in 1000 } 
const int probabilityOfReplacing2ndFittestOf3 = 250; // 1*31 = 31 in 1000 }


int random(int n) 
{
  return (int) (((double) n)*rand()/(RAND_MAX+1.0));
}



void onPeak(int individual)
{
  int distanceToTarget0 = 0; for(int i = 0; i < genotypeLength; i++) if(genotypes[individual][i] != target0[i]) distanceToTarget0++;
  int distanceToTarget1 = 0; for(int i = 0; i < genotypeLength; i++) if(genotypes[individual][i] != target1[i]) distanceToTarget1++;

  // Check whether the current individual is on peak0, and if it is, set peak0 to true
  if(peak0 == false)
  {  
    if(distanceToTarget0 < radius0) 
    {
      peak0 = true;
    }
  }

  // As above but for peak1
  if(peak1 == false)
  {  
    if(distanceToTarget1 < radius1) 
    {
      peak1 = true;
    }
  }
}



void calculateIndividualFitness(int individual)
{
  /* Goes along the length of the current genotype
     Compares each bit with the corresponding position in both the target genotypes
     For every bit that matches the corresponding bit in the target, distanceToTarget is incremented
     Each distanceToTarget will therefore equal the distance away the current genotype is from each target genotype, in terms of bit differences */
  int distanceToTarget0 = 0; for(int i = 0; i < genotypeLength; i++) if(genotypes[individual][i] != target0[i]) distanceToTarget0++;
  int distanceToTarget1 = 0; for(int i = 0; i < genotypeLength; i++) if(genotypes[individual][i] != target1[i]) distanceToTarget1++;


  double f0 = fitness0 * (1 - (distanceToTarget0/(double)radius0));
  double f1 = fitness1 * (1 - (distanceToTarget1/(double)radius1));

  fitness[individual] = 0;

  if(f0 > 0)
  {
    fitness[individual] = f0;
  }

  if(f1 > fitness[individual])
  {
    fitness[individual] = f1;
  }
}





void initialise()
{
  unsigned int s; ifstream u("/dev/urandom",ios::in|ios::binary);
u.read((char*)&s,sizeof(s)); if (!u) {cerr << "Failed to get seed from kernel random number source device" << endl; exit(1);}
  srand(s);

  /* Initialise target0 and target1 (which represent the two peaks)
     target0 stores the genotype at the top of peak 0
     target1 stores the genotype at the top of peak 1 */
  for(int i = 0; i < genotypeLength; i++)
  {
    target0[i] = 0; // target0 is a sequence of 0s
    target1[i] = target0[i]; // target1 is initialised to all 0s but will be mutated according to the value of distanceBetweenTargets in the next step
  }


  /* The following code mutates target1 a number of times equal to distanceBetweenTargets
     This ensures target1 is the correct number of bits different in fitness to target0 */

  /* change stores true or false for each bit, to ensure that once a bit gets mutated it cannot be mutated back
     i.e., preventing back mutation */
  bool change[genotypeLength]; 

  // Initialises each position in the change array to false
  for(int i = 0; i < genotypeLength; i++) 
  {
    change[i] = false;
  }

  /* Changes a number of bits in target1 to make it the required number of bits
     different to target0 */
  for(int d = 0; d < distanceBetweenTargets; d++) 
  {
    int i;
    
    do 
    {
      i = random(genotypeLength); // Randomly chooses a position in the genotype that hasn't already been changed
    } 
    while(change[i] != false);

    change[i] = true; // Sets change to true to ensure this position cannot be selected again
    
    int base;
    
    do
    {
      base = random(alphabetSize); // Randomly generates a new value to go in the selected position of the genotype
    } 
    while(base == target1[i]); // Ensures the new value isn't the same as the one already stored in that position of the target1 genotype

    target1[i] = base; // Replaces the existing value with the mutant
  }



  // Generate the initial population
  for(int i = 0; i < numGenotypes; i++)
  {
    for(int j = 0; j < genotypeLength; j++)
    {
      genotypes[i][j] = target0[j]; // Initialises the whole population to target0
    }

    calculateIndividualFitness(i); // Calculates the initial fitness of each individual in the population
  }

  bool *chosen;
  chosen = new bool[numGenotypes];

  // Initialise chosen to false
  for(int i = 0; i < numGenotypes; i++)
  {
    chosen[i] = false;
  }

  // Puts half of the population on to target1 (so that there is 50% of the population starting on each peak)
  for(int i = 0; i < numGenotypes/2; i++)
  {
    int randomIndividual;

    // Picks an individual at random that hasn't been previously chosen
    do
    {
      randomIndividual = random(numGenotypes);
    }
    while(chosen[randomIndividual]);

    chosen[randomIndividual] = true;

    // Places this individual on peak 1
    for(int j = 0; j < genotypeLength; j++)
    {
      genotypes[randomIndividual][j] = target1[j]; 
    }

    calculateIndividualFitness(randomIndividual); // Calculates the updated fitness of each individual that has been moved to target1  
  }

  delete chosen;
}




/* Copies parent A's genotype up to a certain point, and then copies parent B's genotype
   The resulting sequence replaces the child in the genotypes array */
#warning "Crossover function must be able to cope with child==parentA or child==parentB ... okay for this code for standard crossover"
void crossover(int parentA, int parentB, int child)
{
  // Randomly generates a position in the genotype at which to cut (0 < cutPoint > genotypeLength to ensure crossover occurs)
  int cutPoint = 1 + random(genotypeLength - 1);

  for(int g = 0; g < cutPoint; g++)
  {
    genotypes[child][g] = genotypes[parentA][g]; // Copies parent A's genotype up to the cut point
  }  
  for(int g = cutPoint; g < genotypeLength; g++)
  {
    genotypes[child][g] = genotypes[parentB][g]; // Copies parent B's genotype from the cut point to the end of the sequence
  }
}




/* Mutate each position in the individual sequence with a given probability (mutationRate)
   There is an equal chance of mutation to any of the other possible values (GATC) */
void mutate(int individual)
{
  int r = 0; 
  int newValue;

  // Goes along each position of the sequence in the current individual
  for(int j = 0; j < genotypeLength; j++)
  {
    r = random(1000000);

    // If the mutation rate is greater than or equal to random value r, mutate the current position
    if(mutationRate > r)
    {
      /* Do the point mutation by randomly generating a new value to replace the original value stored in the genotypes array
         N.B. The new value cannot be the same as the original value (otherwise there has been no mutation) */

      newValue = genotypes[individual][j]; // Copies the original value into newValue, which is then mutated

      while(newValue == genotypes[individual][j]) // Makes sure the mutation is done again if necessary
      {
        newValue = random(alphabetSize);
      }

      genotypes[individual][j] = newValue; // Replaces the original value in the genotypes array with the new mutant value
    }
  }
}




void nextStep()
{
  int *popArray;
  popArray = new int[numGenotypes];

  // Initialises popArray as an ordered list of 0 up to numGenotypes-1
  for(int i = 0; i < numGenotypes; i++)
  {
    popArray[i] = i;
  }

  // Sets the initial position to be at the top of the list
  int position = numGenotypes - 1;

  // Goes along the length of popArray and changes the order to random
  for(int i = 0; i < numGenotypes; i++)
  {
    int replace = random(position+1); // Generates a random number between 0 and the current position

    // Switches the value stored in position with the value at the randomly selected position
    int temp = popArray[position];
    popArray[position] = popArray[replace];
    popArray[replace] = temp;

    // Moves position down one in the list
    position--;
  }


  int a,b,c;

  for(int i = 0; i < numGenotypes-2; i+=3)
  {
    a = popArray[i];
    b = popArray[i+1];
    c = popArray[i+2]; 

   // Reorder such that c is the least fit
  if(fitness[a] < fitness[c]) 
  {
    int temp = a;
    a = c;
    c = temp; 
  }
  if(fitness[b] < fitness[c]) 
  { 
    int temp = b;
    b = c;
    c = temp;
  }

  // Make sure a is the fittest
  if(fitness[a] < fitness[b])
  {
    int temp = a;
    a = b;
    b = temp;
  }   

    int child;
#define  nonStandardSteadyStageGA
#ifndef  nonStandardSteadyStageGA
    child = c;
#else  //nonStandardSteadyStageGA
  #warning "Non-standard Steady Stage GA: individual to be replaced chosen by rank based selection where probabilities of selecting fittest and second fittest of 3 may be non-zero"
    int r = random(1000);
    if (r<probabilityOfReplacing1stFittestOf3) {child = a; a = c;} // Copy c into a?
    else if
    (r<(probabilityOfReplacing1stFittestOf3+probabilityOfReplacing2ndFittestOf3)) {child = b; b = c;}
    else child = c;
    // nb may well be that child==a or child==b but standard crossover code (below and in crossover()) can cope with that
#endif

    // Crossover a, b (in random order) to create child that replaces c
    if(random(2) == 1) // Randomly generate either a 2 or a 1. If it's 1, switch a and b.            
    {
      int temp = a; 
      a = b;
      b = temp;
    }
    // else leave a and b as they are

    // Crossover a and b and replace c with the resulting child
    crossover(a,b,child);

    // Mutate child 
    for(int m = 0; m < numMutationsPerReproduction; m++) 
    {
      mutate(child);
    }

    calculateIndividualFitness(child);
  }

  delete popArray;
}




int main(int argc, char* argv[])
{ 
  //int n;

  if(argv[1]==0) {cerr << "Missing argument" << endl; exit(1);}
  numGenotypes = atoi(argv[1]);

  genotypes = new int*[numGenotypes];
  for (int g=0;g<numGenotypes;g++)
     genotypes[g] = new int[genotypeLength];

  fitness = new double[numGenotypes];

  if(argv[2]==0) {cerr << "Missing argument" << endl; exit(1);}
  mutationRate = atoi(argv[2]);

  initialise();

  int peak0LostAt = -1;
  int peak1LostAt = -1;

  for(int reproduction = 1; reproduction <= numGenerations; reproduction++)
  {
    peak0 = peak1 = false;
    nextStep();

    // Works out whether there are any individuals on either peak
    for(int i = 0; i < numGenotypes; i++)
    {
      onPeak(i);
    }
   
    // If it's the first time the peak has been lost, record the current reproduction
    if ((peak0==false) && (peak0LostAt==-1)) peak0LostAt = reproduction;
    if ((peak1==false) && (peak1LostAt==-1)) peak1LostAt = reproduction; 
  
    // If both peaks have been lost, output when they were lost and end the program
    if((peak0LostAt > -1) && (peak1LostAt > -1)) break;
  }

  cout << "Population size = " << numGenotypes << ", mutation rate =  " << mutationRate << ", peak 0 (fit) peak lost at =  " << peak0LostAt << ", peak 1 (flat) peak lost at =  " << peak1LostAt << endl;

  delete fitness;
  for (int g=0;g<numGenotypes;g++)
     delete genotypes[g];
  delete genotypes;
  
  return 0;
}

